/******************************************************************************************
  This is to plot the moments and profiles from baseline model and other robustness check models
  in the same graph.

  in their own folder, 

  ** dta_MSM_Moments.dta: moments
        _f: generate from Fortram program
        _d: from SIPP data
       lfpr lnw sd_lnw lnw_fd lfpr_diff
*******************************************************************************************/

*by Xiaodong Fan, fanxiaodong@gmail.com

#delimit ;

cap log close;
clear all;
drop _all;
set more 1;
pause on;

global gpicsdir "pics_rob_compare";
global gdtadir ".";

include do_globals.do;

log using ${gpicsdir}/log_MSM_rob_compare.log, replace;

*********************************************************************************************************************;
*  Moments from baseline and various models;
*  1. consumption floor larger-say 2.5
*  2. consumption floor lower say 1.8
*  3. I would move discount rate to 0.95 and change r to 0.05 at same time
*  4. we ought to see what would happen if consumption were to decline with age. Keep r at 1.03 and decrease beta to 0.95 to see what happens
*  5. If we do initial wealth it should be a smaller number, say 50,000 which still seems like a large amount of money for young guys
*  6. C and SSA has no age trend, not matching C or SSA moments either;
*********************************************************************************************************************;
local lvdir0 ".";
local lvdir1 "../../../uBP_MPI_31/main/output";
local lvdir2 "../../../uBP_MPI_32/main/output";
local lvdir3 "../../../uBP_MPI_33/main/output";
local lvdir4 "../../../uBP_MPI_34/main/output";
local lvdir5 "../../../uBP_MPI_35/main/output";
local lvdir6 "../../../uBP_MPI_noage_noms/main/output";

*local lvrobnums 1 2 3 4 5;
local lvrobnums 1 2 3 5 6;

*************************************************************************;
insheet using pars_mapping.csv, comma names clear;
sort incode;
save dta_temp.dta, replace;

* get parameter estimates;
foreach iv in 0 `lvrobnums' {;
   infile par using `lvdir`iv''/txt_parin9.txt, clear;
   if (`iv'==6) {;
       keep if [_n]<=19;
       gen incode = [_n]; 
       replace incode = incode + 3 if incode >= 17;
       replace incode = incode + 3 if incode >= 23;
   };
   else {; 
       keep if [_n]<=${gvparnum};
       gen incode = [_n];
   };
   destring par, replace;
   prog_parSSAx1000 par;

   gen byte model = `iv';
   if (`iv' > 0) {;
      append using ${gdtadir}/dta_rob_pars.dta;
   };
   sort incode model;
   save ${gdtadir}/dta_rob_pars.dta, replace;
};
merge m:1 incode using dta_temp.dta;
keep if _merge==3;
drop _merge;

sort inpaper model;

reshape wide par, i(inpaper) j(model);
sort inpaper; 

local lvpars par0 par1 par2 par3 par5 par6;
*prog_parSSAx1000 `lvpars';
export excel inpaper description `lvpars' incode using "table4paper", sheet("pars_rob") sheetreplace firstrow(variables);

*************************************************************************;
* import moments for comparison;
* from the robustness check models only, not including the baseline model;
foreach iv of local lvrobnums {;
   if (10<1) {;
      prog_MSM_moments_infile "`lvdir`iv''/txt_parin9_mindis_moments.txt" iter;
      sum iter;
      local lviterm = r(max);
      keep if iter==`lviterm';
      drop iter;
   };
   else {;
      prog_MSM_moments_infile "`lvdir`iv''/txt_MSM_Moments.txt" "";
   };
   ren *_f *; 
   sort t;

   gen byte model = `iv';
   if (`iv' > 1) {;
      append using ${gdtadir}/dta_MSM_Moments_rob_compare.dta;
   };
   sort model t;
   save ${gdtadir}/dta_MSM_Moments_rob_compare.dta, replace;
};


local lvmoments lfpr lnw sdlnw lnw_fd alfpr1to0 alfpr0to1 C ssa plnw41;

keep t model `lvmoments';
reshape wide `lvmoments', i(t) j(model);
sort t;
save ${gdtadir}/dta_MSM_Moments_rob_compare.dta, replace;


use ${gdtadir}/dta_sipp_Moments.dta, clear;
sort t;
merge 1:1 t using ${gdtadir}/dta_MSM_Moments_rob_compare.dta, nogenerate;
sort t;
save ${gdtadir}/dta_MSM_Moments_rob_compare.dta, replace;

sort t;

* plot graph;

foreach iv in lfpr lnw sdlnw lnw_fd C ssa {;
   local lyange "";
   local lylabel "";
   if ("`iv'"=="lfpr_diff") {;
      local lvtperiod if t>=${gvdata_diff} & t<=${gvdata10};
      local lvxlabel "xlabel(${gvdata_diff}(5)60 62 ${gvdata10},grid)";
      local lvytitle "Difference in LFPR";
      local lvpos = 6;
   };
   else if ("`iv'"=="ssa") {;
      local lvtperiod if t>=${gvdata_ss} & t<=${gvdata10_ss};
      local lvxlabel "xlabel(${gvdata_ss}(1)${gvdata10_ss},grid)";
      local lvytitle "Social Security Application";
      local lvpos = 1;
   };

   else {;
      local lvtperiod "if t>=${gvdata0} & t<=${gvdata10}";
      local lvxlabel "xlabel(${gvdata0} 30(10)60 62 ${gvdata10},grid)";
      if ("`iv'"=="lfpr") {;
         local lvytitle "Labor Force Participation Rates";
         local lvpos = 7;
      };
      else if ("`iv'"=="lnw") {;
         local lvytitle "Log Wages";
         local lvpos = 6;
      };
      else if ("`iv'"=="lnw_fd") {;
         local lvytitle "Log Wages (FD)";
         local lvpos = 6;
      };
      else if ("`iv'"=="sdlnw") {;
         local lvytitle "Log Wages: Standard Deviation";
         local lvpos = 11;
         local lyange 0 1;
         local lylabel 0(0.2)1;
      };
      else if ("`iv'"=="C") {;
         local lvytitle "Consumption";
         local lvpos = 6;
         local lyange 4 12;
         local lylabel 4(2)12;
      };
   };


   line `iv'_d `iv'1 `iv'2 `iv'3 `iv'5 `iv'6 t `lvtperiod', lpattern("l" "_" "_." "-" "-." "__--" "_.-." ) graphregion(color(white))
        `lvxlabel' ytitle("`lvytitle'") yscale(range(`lyrange')) ylabel(`lylabel')
        legend(pos(`lvpos') ring(0) col(1) order(1 "Data" 2 "Larger CF" 3 "Smaller CF" 4 "Larger r" 5 "Larger A0" 6 "No age trend"));
   graph export ${gpicsdir}/eps_moments_rob_compare_`iv'.eps,replace;


};


cap log close;



